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Abstract 


In  this  first  of  two  papers  ,  computable  a  posteriori 
estimates  of  the  space  discretization  error  in  the  finite 
element  method  of  lines  solution  of  parabolic  equations  are 
analyzed  for  time-independent  space  meshes.  The  effectivity 
of  the  error  estimator  is  related  to  conditions  on  the  solu¬ 
tion  regularity,  mesh  family  type,  and  asymptotic  range  for 
the  mesh  size.  For  clarity  the  results  are  limited  to  a 
model  problem  in  which  piecewise  linear  elements  in  one 
space  dimension  are  used.  The  results  extend  straight¬ 
forwardly  to  systems  of  equations  and  higher  order  elements 
in  one  space  dimension,  while  the  higher  dimensional  case 
requires  additional  considerations.  The  theory  presented 
here  provides  the  basis  for  the  analysis  and  adaptive  con¬ 
struction  of  time-dependent  space  meshes,  which  is  the  sub¬ 
ject  of  the  second  paper.  Computational  results  show  that 
the  approach  is  practically  very  effective  and  suggest  that 
it  can  be  used  for  solving  more  general  problems. 


\ 


1.  Introduction 


In  recent  years  interest  has  grown  in  the  method  of  lines 
(MOL)  approach  for  the  numerical  solution  of  time-dependent 
partial  differential  equations  (PDEs)  arising  in  biology,  rhe¬ 
ology,  structural  and  fluid  mechanics,  and  many  other  fields. 
In  this  approach  a  problem  is  first  discretized  in  space,  for 
example,  by  a  Galerkin  or  difference  method.  This  results  in 
a  system  of  ordinary  differential  equations  (ODEs)  which  can 
be  efficiently  solved  by  one  of  the  available  variable-step, 
variable-order  ODE  packages.  This  segmented  approach  is  obvi¬ 
ously  practical,  and  modern  ODE  integrators  are  generally  re¬ 
liable  in  estimating  and  controlling  the  error  due  to  the  time 
discretization.  However,  one  is  usually  interested  in  the 
total  error,  i.e.  the  error  of  the  approximate  solution  with 
respect  to  the  exact  solution  of  the  original  PDE.  The  advan¬ 
tages  of  highly  accurate  ODE  integrators  are  diminished  if  the 
error  due  to  the  space  discretization  is  large.  Those  avail¬ 
able  programs  implementing  both  segments  of  the  MOL  approach 
(cf.  survey  in  [16])  generally  have  no  facility  to  estimate 
and  control  both  components  of  the  error. 

Recently  a  MOL  procedure,  which  we  shall  call  the  finite 
element  method  of  lines  (FEMOL),  was  proposed  for  the  solution 
of  parabolic  PDEs  (cf.  [1]).  Piecewise  linear  finite  elements 
in  one  space  dimension  were  employed  in  a  model  problem,  which 
consisted  of  a  second  order  homogeneous  PDE,  with  smooth  and 
compatible  initial  data.  Using  the  theory  and  practice  of  a 


posteriori  error  estimates  developed  for  elliptic  PDEs  in 
[2],  [3],  [4],  and  [5],  an  estimator  for  the  space  discretiza¬ 
tion  error  was  constructed.  This  estimator  and  a  novel  adap¬ 
tive  time  discretization  scheme  were  utilized  to  estimate  the 
total  error  in  the  FEMOL  in  the  case  of  time-independent  space 
meshes . 

We  shall  restrict  our  attention  here  to  the  space  discre¬ 
tization  error  in  the  FEMOL,  i.e.  the  total  error,  assuming 
that  the  resulting  ODE  system  can  be  solved  exactly.  In  any 
given  application  this  of  course  cannot  be  accomplished,  and 
the  effects  of  the  errors  in  the  ODEs  must  be  taken  into 
account.  The  purpose  of  this  paper  is  to  completly  isolate  the 
space  component  of  the  error  and  means  of  estimating  it.  By 
keeping  the  analysis  independent  of  the  time  discretization 
scheme,  none  of  the  state  of  the  art  ODE  solvers  are  excluded 
a  priori  when  dealing  with  the  important  question  of  how  to 
properly  balance  the  two  components  of  the  error. 

In  this  first  of  two  papers,  it  is  assumed  that  the  space 
mesh  remains  fixed  throughout  the  time  evolution  of  the  problem. 
Using  piecewise  linear  elements  in  one  space  dimension we  consider 
the  model  problem  of  a  second  order  nonhomogeneous  parabolic  PDE  ,  with 
initial  data  not  necessarily  smooth  or  compatible  .  For  related  results 
with  higher  order  elements  in  the  setting  of  coupled  systems  of 
equations  in  one  space  dimension  see  [8]. 

The  results  given  here  are  extended  in  a  second  paper  [7] 
to  analyze  means  of  estimating  the  space  discretization  error 


when  the  space  mesh  is  allowed  to  change  discontinuously  in 
time.  Based  upon  these  estimates,  a  procedure  is  given  there 
for  the  control  of  the  space  discretizaticn  error  by  the  adap¬ 
tive  construction  of  space  meshes  during  the  solution  process. 
This  algorithm  is  then  tested  in  some  computational  examples. 

The  effectivity  of  the  space  discretization  error  estima¬ 
tor  constructed  in  [1]  depends  on  properties  of  the  exact  solu¬ 
tion  in  the  space-time  domain  and  is  asymptotic  in  nature.  We 
give  sufficient  conditions  and  a  priori  bounds  on  the  asympto¬ 
tic  range  for  this  estimator  to  work.  These  bounds  are  overly 
pessimistic  for  any  given  problem,  and  we  shall  see  in  some 
examples  that  the  results  are  much  better  than  the  theory  pre¬ 
dicts.  Under  less  restrictive  conditions  and  for  all  mesh 
partitions,  it  is  also  shown  that  a  modification  of  our  prin¬ 
cipal  error  estimator  can  be  expected  to  perform  well.  This 
modification  takes  more  fully  into  account  the  case  where  an 
exact  solution  has  a  strong  time-dependence  with  respect  to 
the  size  of  the  mesh  used. 

This  paper  consists  of  the  following  sections.  In  section 
2  notation,  mathematical  preliminaries,  and  the  model  problem 
are  introduced.  In  section  3  certain  regularity  classes  are 
defined  and  an  accompanying  theorem  is  proven,  showing  that 
the  model  problem  is  well-posed.  Based  upon  this  proof  and 
result,  in  section  4  a  priori  estimates  needed  in  section  5 
are  given.  Our  a  posteriori  error  estimates  are  presented  in 
sections  5  and  6.  Section  7  consists  of  the  description  of 
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computational  procedures  used  and  examples. 

Many  of  the  required  function  spaces,  results,  and  tech¬ 
niques  used  in  sections  2-6  to  prove  our  results  are  somewhat 
technical.  In  order  to  facilitate  a  first  reading  of  the 
paper,  we  present  here  some  basic  notions  and  a  sample  of  our 
computational  results.  By  then  skipping  directly  to  section 
7,  the  reader  should  have  a  good  overview  of  the  main  ideas 
presented . 

Let  u  =  u(t,x)  be  the  solution  of  the  model  problem 


( u  =  f-Lu  =  f(t,x)+[a(x)u  ]  -b(x)u; 

IT  XX 

Cl . 1)  (  u(t ,x)  =  g(t  ,x) ; 
ku( o ,x)  =  uQ(x) ; 


0  <  x  <  1 ,  t>0, 

x  =  0 ,1 ,  t  >  0  , 
0  <  x  <  1 , 


where  a>0,  b  >  0,  uQ,  g,  and  f  are  given 
functions . 


The  FEMOL  solution  of  eqs.  (1.1)  is  accomplished  by  parti¬ 
tioning  the  interval  (0,1)  according  to 

(N)  .  . 

A  =  {0  =  xQ  <  x^  <  • • •  <  x^  =  1}  and  determining  the  unique 

function  U  =  U(t,x),  which  is  piecewise  linear  in  x  and 

smooth  in  t,  and  satisfies  a  weak,  or  integral  form  of  eqs. 

(1.1).  For  each  t,  U  is  characterized  by  the  vector  of  its 

values  {U(t,x.)>.  ,  ..  ,  ,  which  is  the  solution  of  an  N-l 

dimensional  ODE  system,  assumed  here  to  be  exactly  solvable. 

In  the  computational  examples,  the  ODE  systems  were  solved  by 

using  successively  smaller  time  discretization  error  tolerances, 


5 


the  process  being  terminated  when  all  relevant  quantities 
experienced  only  insignificant  changes  in  the  time  intervals 
of  interest. 

For  each  t  >  0  the  goal  is  to  estimate 
(1.2)  HI  e  ( t ,  -  )(||  =  f  a(x)e2(t,x)dx 

U  o  x 

where  e  =  u  -  U  is  the  space  discretization 
error . 

In  order  to  estimate  |||e(t,*)|||  ,  we  define  the  local  error 

indicators  {^.(t)}.  ,  ..  by 
J  J  ■ 1 


(1.5) 


0(t) 


F(t)/|||e(t,  •  )!l!  ; 


t  >  0  . 


r 


£(•)  is  an  effective  estimator  if  there  exists  a  reason¬ 
able  constant  C  (depending  on  the  admissible  solution  class 

for  eqs.  (1.1)  and  the  admissible  class  of  mesh  partitions  to 
( N) 

which  A  belongs,  but  not  on  the  magnitude  of  the  data 

uQ  ,  f,  or  g)  such  that 


(1.6) 


l/0(t)  <  C; 


Vt  >  0. 


Moreover,  it  is  desirable  that 


(1.7)  0(t)  -+■  1  as  the  partitions  of  (0,1)  are  refined. 


In  one  of  the  examples  of  section  7 ,  the  exact  solution  of 
eqs.  (1.1)  is  taken  to  be 


(1.8) 


u(t,x)  =  |  |  tanh  [2p(x-10t) ] ;  p  =  20, 


which  is  a  wave  with  approximate  front  width  p  that  moves 
ir.  the  positive  x-direction  at  speed  10  (cf.  figure  1(a), 
(b),(c)).  For  each  of  various  values  of  N,  the  FEMOL  was 
implemented  in  this  example  by  defining  the  partition  A  ^  ^ 
according  to  Xj  =  j/N;  j  =  0,N.  The  tables  below  illustrate 
for  various  t  the  performance  of  E(t)  as  compared  with  the 
relative  error  EREL(t),  defined  by 


(1.9) 


EREL(t) 


TABLE  2.  0(t) 


20 

-  - 

40 

80 

-  ; 

160  I 

1.320 

.995 

.989 

.996 

.776 

1.076 

.989 

.996 

.806 

.903 

.988 

.996 

1.685 

.961 

.987 

.996 
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X  X  X 

(a)  (b)  (c) 


Figure  1.  Solution  u(T,X)  in  Eq.  (1.8)  vs.  X  for  various  T. 

The  last  three  columns  of  Table  1  illustrate  for  each  t 
the  expected  linear  rate  of  convergence  for  |||  e  ( t ,  •  )|||  .  For 
N  =  20  the  mesh  used  is  not  fine  enough  to  resolve  the  wave- 
front  and  the  errors  oscillate  about  the  50%  mark. 

The  numbers  in  Table  2  clearly  illustrate  the  effectivity 
of  estimating  |||e(t,*)|||  with  E(t).  While  values  for  the 
constant  C  in  (1.6)  on  the  order  of  10  or  more  would  be 
acceptable  in  many  applications,  for  none  of  the  entries  in 
Table  2  is  such  a  pessimistic  bound  realized.  We  see  that 
for  about  13%  or  less  error,  the  estimator  £(•)  differs 
from  HI e (  -  ,  *  )|||  by  no  more  than  about  1%,  a  truly  remarkable 
result . 

These  and  other  computational  results  are  discussed  more 
fully  in  section  7.  We  again  mention  that  the  estimator 
£(•)  can  be  used  in  an  algorithm,  where  decisions  concerning 
modifications  of  the  space  mesh  are  made  adaptively  during  the 
solution  process.  Provided  that  proper  constraints  are  imposed, 


this  then  leads  to  an  implementation  of  the  FEMOL  which  is 
robust  enough  to  adaptively  modify  meshes  in  the  case  of  local¬ 
ized  activity,  such  as  for  the  travelling  wave  above,  and  yet 
is  still  efficient  enough  to  handle  problems  where  little  or 
no  mesh  modification  is  necessary. 
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2 .  Mathematical  Preliminaries  and  the  Model  Problem 

Throughout  this  paper  we  shall  use  the  convention  that 
the  variables  i,  j,  l ,  m,  N,  p,  and  q  take  only 
integer  values. 

Let  J  be  a  bounded  open  interval  in  IR^  ,  J  its  clo¬ 
sure,  and  V  be  an  arbitrary  Hilbert  space.  We  denote  by 
C°(J;V)  those  V-valued  functions  which  are  continuous  on  J, 
C°°(J;V)  the  subset  of  infinitely  differentiable  functions  on 
J  for  which  all  derivatives  have  continuous  extensions  on  J, 
and  Cg(J;V)  those  functions  in  C“(J;V)  with  compact  sup¬ 
port  in  J.  When  V  =  IR^  we  suppress  V  and  write  C°(J), 
C"(J),  C"(J). 

Let  I  =  (0,1)  be  the  unit  open  interval  in  3R1,  I  its 
closure,  and  [0,T)  be  a  bounded  half  open  interval  in  IR'*'  . 
For  4>  €  C  ( I )  ,  set 


«t> 


r 

$(x)ty(x)dx  and 
I 


<  <f>  ,4» . 


For  each  i  2  0  the  Sobolev  space  H1{Hq}  is  the  completion 
of  C°°(I)  {Cq(I)  }  with  respect  to  the  norm 


Ml  =  l 

p  =  0 

and  for  each  real  s  >  0,  Hs  with  norm  ||  •  ||  is  a  Hilbert 

b 

space  which  is  defined  by  interpolation  as  in  [15]. 

Let  a  and  b  be  sufficiently  smooth  functions  on  I 


for  which 
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(2.1) 


(2.2) 


0  <  a  5  a(x)  5  a 


0  s  b(x)  <  b 


Vx  6  I. 


For  u,v  €  Hn  we  define  the  bilinear  form 


(2.3) 


B(u,v)  =  <au  ,v  >  +  <bu,v>. 

’  X  X  ’ 


By  (2.1)  and  (2.2)  3  positive  constants  C^,C2  suc^ 


that 


(2.4) 


C2||w||^  <  B(w,w)  £  C1||w!|  J 


Vw  €  Hj, 


and  i| w If E  =  B(w,w)  defines  a  norm  equivalent  to 


The  operator  L  defined  by 


Lu  =  -(au  )  +  bu 

X  X 


2  1 

with  domain  V(L )  =  H  fl  is  the  unique  positive  self- 
adjoint  operator  induced  by  B.  Using  spectral  representations 
(cf.  [9],  [13],  [14])  we  can  define  fractional  powers  of  L. 
Results  in  [11]  show  that  for  a  2  0,  a/2  *  integer  +i 


(2.5) 


D(La/2)  =  {v  €  Ha  :  LPv(x) |x=0il  =  0; 


V p  satisfying  0  £  p  <  a/2  -  1/4} 


and,  as  in  (2.4), 


(2.6)  the  spectral  norm  ||L 


Q  is  equivalent  to  IHIa 

nr-!  a/2\ 

on  P(L  ) . 
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Let  P 

be  a  family  o 

f  mesh 

partitions  of  I  and 

A 

=  {0  =  xQ  <  x 

1  < 

<  xn(a)  =  6  p‘ 

For  j 

=  1 ,  N  (  A  )  we 

write 

I. 

=  (x .  , ,x . ) , 

h. 

=  x .  -  x .  , ,  and 

2 

3-1  3 

3 

3  3-1 

h(  A ) 

=  min  h .  , 

h(A)  =  max  h. 

3 

lsjsN( A) 

3 

l<j<N( a) 

P  is 

said  to 

be  a 

«(X)-regular  family  if  3  constants 

<  >  1 

o 

A 

such  that 

(2.7) 

h( A )  >  Xh^CA) 

VA 

€  P. 

For  v 

w  €  Hj, 

j  =  1 , N ( A ) ,  and  +  €  C°(T.) 

we 

shall 

write 

*  J 

fx. 

(x)dx , 

3 

► 

A  A 

j 

•j 

d 

J  J 

Ca(x)wx(x)  +  b(x)w  (x)]dx,  and 

For  each  A  €  P, 

whose  restrictions  to 
A  A 

PQ»  pl5  mapping 

S(A)  are  defined  by 


S  ( A )  c  H 


any 


Q  is  the  subset  of  functions 
are  linear.  The  operators 
S ( A ) ,  respectively,  to 
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'  <P*v-v,4»  =  0 

(2.8)  \  B( P^v-v ,  4> )  =  0  V*  €  S(A). 

v  <L^v,<t»  =  B(.v  ,<t> ) 

is  a  positive  self-adjoint  operator  on  S(A)  satisfy¬ 
ing 

(2.9)  lA<P1u)  =  Po(Lu)  Vu  €  P(L)* 

By  (2.1)  and  (2.2)  3  a  positive  constant  C  such  that 

V<f>  €  S  (  A ) 

(2.10)  5  CllL”1*  |IE  5  c2IU  ll0, 

and  if  p  is  k  ( M -regular 

(2.11)  II  la^  II 0  '  ( C/^)h-K  (  A  )!)<#>!!  E . 

The  following  lemma  is  well-known  (cf.  [6,  Chap.  4]). 

Lemma  2.1.  There  exists  a  positive  constant  C  such 
that  for  each  A  €  p 

(2.12)  ll(l-pj)vlls  S  Cha_S(A)||vlla;  Vv  €  P(La/2), 

0Ss-1Sct£2j 

(2.13)  II  (  I-Pq  )  vll  Q  2  Cha(A)llvllCT;  Vv  €P(LCT/2), 

0  <  a  <  2,  °  *  1/2 . 


We  shall  consider  the  following  model  problem. 
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^ ~  u(t,x)  +  Lu(t,x) 


(2.14)  \ 


u( t , 0 )  =  0 

u(0 ,x) 


f(t,x) ;  t,x  €  (OjT)  *  I, 
u(t  ,1) ;  t  €  ( 0 ,T) , 
uQ (x) ;  x  €  I . 


The  data  and  f  are  such  that  the  solution  . 

u(;uQ,f)  of  eqs.  (2.14)  resides  in  a  certain  regularity  class 
H,  which  will  be  specified  in  section  3. 

The  weak  form  of  eqs.  (2.14)  is:  find  u  :  (0,T)  -+•  Hq 
satisfying 


'<u  (t)  ,<#»>+  B(u(t)  ,*)  =  <f  (t)  ,<*»>;  Vt  ,*  €  (  0  ,T)  x  hJ, 


(2.15) 


u(0)  =  Uq • 


For  each  A  ?  P,  the  FEMOL  solution  U  :  [0,T)  ■*  S(A) 
is  defined  by 


(2.16) 


'<U*(t)  ,<*»  +  B(UA(t)  ,*)  =  <f(t),*>;  €  ( 0  ,T)  x  S(A), 

UA(0)  =  uj  €  S ( A ) , 


where  is  chosen  to  approximate  u^  and  is  specified  in 

section  4.  Eqs.  (2.16)  may  be  rewritten  as 


(2.17) 


’UA(t)  +  LAUA(t)  =  Pgf(t); 


t  €  (0,T)  , 


.lAo) 


U0’ 


and  for  reasonable  data  u^ ,  f  are  equivalent  to  a  uniquely 
solvable  N(A)  -  1  dimensional  ODE  initial  value  problem. 
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The  following  discrete  smoothing  property  can  easily  be 
verified  (cf.  [9]). 

Lemma  2.2.  Let  =  U^(;Uq,0)  denote  the  solution  of 

eqs.  (2.17)  with  P^f  =  0.  Then  3  a  positive  constant  C 
such  that  for  i  =  0  or  1 


(2.18) 


A 

U  (t) 


<  Ct-(i+2j)/2 


liu£|| 


0"0 


Vt  >  0, 


>  0. 


Remark.  When  there  is  no  threat  of  confusion,  we  shall 
often  drop  a  when  writing  S,  h,  PQ ,  P1  ,  U,  and  UQ . 

Our  work  requires  the  use  of  the  nonisotropic  Sobolev 
2  s  s 

spaces  H  ’  ,  which  we  now  define.  For  nonnegative  i  and 
m,  let  H^(J;Hm)  be  the  completion  of  C°°(J;Hm)  with  res¬ 
pect  to  the  norm 


J 

H*(Hm) 


dt , 


and  define  H2*’*(J)  =  H°(J;H2*)  n  H*(J;H°)  with  norm 

2  s  s  J 

For  each  real  s  >  0  the  space  H  ’  (J)  with  norm  ||*||_ 

/  S  j  5 

is  a  Hilbert  space  which  is  defined  by  interpolation  as  in 
[15],  We  state  in  convenient  forms  some  well-known  trace  and 
embedding  results. 


1 
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Theorem  2.1.  (cf.  [15,  Thms .  1.2.3,  1.4.1,  Prop.  4.2.3]) 
Let  s  >  0,  v  €  H2s,S(J),  and  j  satisfy  0  £  j  <  s. 

Then  v  €  (s- j  ) ,  (s- j  )  ^  j ^  and  3  a  positive  constant  C 

such  that 


J 

2 ( s- j )  ,  ( s  — j ) 


s  C  ||  v  || 


J 

2s  ,s 


Theorem  2.2.  (cf.  [15,  Thms.  1.3.1,  1.4.2]) 

Let  s  >  1/2,  v  €  H2s,S(J),  and  j  €  [0,s-l/2).  Then 
v  €  C°(J;H2(s_j'1/2))  and  3  a  positive  constant  C  such 
that 


sup 

t€J 


ftj  V(t 


2 ( s- j -1/ 2  ) 


< 


C  ||  v  || 


J 

2s  ,s 


Theorem  2.3. 

(cf.  [15, 

,  Thm. 

4.2.1,  Prop.  4.2, 

.2]) 

Let  s  >  3/4, 

v  €  H2s’ 

,S(J), 

t*  (  J,  and  j 

and  p 

satisfy  0  <  j  +  p/2  <  s  -  3/4 .  Then 


v(t ,x) 


x=0 ,1 


t  =  t* 


ap  /a 


ax^atJ 


^jV(t.x) 


ft.) 


x=0 ,1 


We  conclude  this  section  by  noting  a  useful  partial  ex¬ 
tension  of  Lemma  2.1. 


Lemma  2.3.  There  exists  a  positive  constant  C  such  that 
for  each  &  (  P 


<  Chau)|!v|| 


J 

2m- 


9 


u  =  i 


v v  €  h2m-’^(j)  n 

+  a/2 ,  j  >  0  , 


H0(J;Hj), 

1  2  a  5  2. 
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Proof.  Set  V  =  H°(  J\V(La/2  )  ) .  Since  c“(  J-,V(Lo/2  )  ) 
is  dense  in  the  Hilbert  space  V ,  I  -  P-^  can  be  continuously 
extended  to  an  operator  mapping  V  to  H°,0(J).  By  (2.12)  we 
have 


(2.19) 


(I-Pj)w 


|J 

*0,0 


<  Cha(A)|iw| 


h°(P(lct/2)) 5 


Vw  €  I ! . 


It  follows  from  the  assumptions  of  the  lemma,  (2.5),  and  Thms . 

a3 

2. 1-2.  3  that  w  =  — j  v  €  l/. 

The  application  of  (2.6)  and  Thm.  2.1  in  (2.19)  completes 
the  proof. 


*  t  ute iA&S  ~  ii 
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3 .  Parabolic  Regularity 

Throughout  the  remainder  of  this  paper,  a  and  j3  denote 
arbitrary  but  fixed  real  numbers  related  to  the  data  in  eqs . 
(2.14)  and  satisfy 

(3.1)  a  >  1,  -1  5  p  £  a,  -3,  a,  a/2  t  integer  +  1/2. 


For  real  a,  [ a]  denotes  the  integral  part  of  a  and 

T'S-  =0  for  q  <  p. 

-  ]=P 

We  begin  by  defining  for  v  €  Hh+1,  g  €  H a’  ^  (  0  ,T ) 
the  p-order  compatibility  measure 


(3.2) 


(v,g) 


[(3/2  +  1/43  , 

I  V  |L  v(x) 

x=  0 , 1  l-Q 


+  f  ’  (~1)P  L*‘P 

P  =  1 


p-i 

atf4g(0’x) 


5 


which  by  Thms .  2.1  -  2.3  is  well-defined. 

The  following  facts  concerning  solutions  u(;tp,f)  of 
eqs.  (2.14)  with  initial  data  $  and  right  hand  side  :  are 
well-known . 


Theorem  3.1.  (Existence;  cf.  [15,  Thm.  3.4.1]) 

Let  *  €  H°,  f  6  H0,0(0,T).  Then  eqs.  (2.14)  have  a 
unique  solution  u(;<p,f)  €  H^((  0  ,T)  ;Hq  )  . 

Theorem  3.2.  (Regularity;  cf.  [15,  Thms.  4.53  ,  4.62]) 

Let  <p  €  Ha+1,  f  €  Ha,a/2(0,T),  and  M^RU,f)  =  0. 

a  +  2,  (a  +  2)/2(Q^T)  and  3  a  positive 


Then  u( ; y, f )  €  H 
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constant  C  such  that 


Hull 


( 0  ,T ) 

a+  2 , ( a+  2 ) / 2 


-  C{!^'!a+1  + 


Ilf1! 


( 0  ,T) 
a  ,  a/  2 


Lemma  3.1.  (cf.  [13,  p.  489],  [14],  and  (2.5)) 

Let  f  s  0  and  «  €  P( L( ^  +  1 } 7 2 } )  .  Then  Vt  >  0 
u(t;p,0)  €  V(L^a  ^/'2)  and  3  a  positive  constant  C  such 
that 

H  u<  t  ;tp  ,  0 )  I!  a  +  1  <  C  t“(a  +  1)/2!l4>!|p  +  1. 

We  therefore  have 


Theorem  3.3.  (Parabolic  Smoothing) 

Let  o,  u(;ip50)  be  as  in  Lemma  3.1  and  0  <  T.,{  <  T. 

Then 

(3.3)  u(;ip,0)  €  Ha  2  ’  ^ a+  2  ^  ^  2  (Tft  ,T)  and  3  a  positive 

constant  C  such  that 


(3.4) 


II 


Proof.  By  Lemma  3.1  u(Tftjq>,0)  €  V(  L(  a+1 }  '  2 )  ,  which 
CR 

is  equivalent  to  ( u(T.,.  ;<p  ,  0  )  ,0)  =  0.  By  applying  Thm.  3.2 

to  the  interval  (T,,{,T)  we  have  (3.3)  and 

(T„,T) 

||u(  ;u(Tft;<p,0)  ,0)  II  a+  2 ,  ( a+  2 )  /  2  5  C''U^T* ’,(p  ’  '^a  +  1 

by  Lemma  3.1  s  C  Tl (  a  +  1 )  '  2 ||<p  ||p  +  1 . 
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We  now  define  certain  regularity  classes  for  data  and 
solutions  of  eqs .  (2.14)  in  order  to  combine  Thms .  3.2,  3.3 


for 

use  in 

Section 

4 .  Set 

(3. 

5) 

5  (0 

T1  \ 

5  ±  > 

.  ,  „P  +  1  ua , a/ 2 

=  { v  ,  g  €  Hr  x  H  ’ 

(  0  ,  T )  :  M 

CR, 

P  (V’ 

,g)  =  0}. 

For 

V 

,g 

€  ID2 

,a  ( 0 ,T )  and  0  £  T 

...  <  T 

set 

(3. 

6) 

a 

(a,T, 

.>. » T  ,  v  ,  g )  =  l|vfla  +  1 

+  «««a° 

jT) 

a/2* 

and 

for 

V, 

g  e 

ma’ 

^  ( 0  ,T)  with  p  <  a 

and 

0  <  Tj,  <  T  set 

(3  . 

7)  C 

a 

(p 

,T* 

»T,v 

,g)  =  (1 + T*"(a+1)/2 

)(MCR(0 

a 

»g>  + 

11  vilp  +  l) 

+  »*C/2- 

For 

0 

< 

T*  < 

T  define 

(3. 

8)  Ha 

,P 

(T* 

,T) 

=  {u  €  Ha+2,<:a  +  2)/2(T 

... ,  T )  :  u 

=  u(  : 

;v,g) 

solves  (2.14)  for 

some 

v,g  < 

:  JDa,?  (  0  ,T)  }  , 

and 

for 

u  ( 

;v, 

g)  € 

Ha,i3(T4,T)  set 

(3. 

9) 

ca(p 

,Ts,T,u)  =  7a<P,T4 

>T,v,g) 

. 

By  properly  decomposing  given  initial  data  we  shall  show 

Theorem  3.4.  Let  0  5  T4  <  T  and  Ug,f  ^  Da’^  (0,T), 
where  p  =  a  if  T*  =  0 .  Then  eqs.  (2.14)  have  a  unique 
solution  u(  ;  Ug  ,  f )  €  Ha’^(T,..,T)  and  3  a  positive  constant 


C  such  that 


21 


(3.10)  il  u: 


(T*,T) 

a+2,(a+2)/2 


5  C  Ca(p,T*,T,u  ,f)  =  C  Ca(p,T...,T,u)  . 


Remark .  The  decomposition  used  in  the  proof  of  Thm. 
3.4  will  be  used  to  demonstrate  the  a  priori  estimates  of 
Section  4. 


Proof .  The  outline  of  the  proof  is  as  follows.  Through 

a  sequence  of  elliptic  boundary  value  problems  we  construct 
1  2 

functions  uQ  and  u^  with  the  following  properties: 


(3.11) 


1  *  2 

U0  U0 


u 


O’ 


(3.12) 


1  ua  +  1  J  ,/CR,  1  _  n 

uQ  €  H  and  (uQ,f)  -  0, 


(3.13) 


2  ,  nn  (p+l)/2. 
uQ  €  P(L  ) , 


and  for  which  3  positive  constants  C, ,  such  that 


0.14)  :iu0!iatl 


(3.15)  llu0!|?tl  <  C2 


■'  U0  ^  a  +  1  ’ 

if 

P  = 

a  , 

M^R(0,f), 

if 

P  < 

a  , 

0, 

if 

(3  =  a , 

MCR( 0 , f )  + 
a 

■luo' 

'p  +  1* 

if 

j3  <  a . 

By  linearity  the  solution  u(;uQ,f)  (whose  existence  is 
given  by  Thm.  3.1)  can  be  decomposed  as 


A 
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(3.16)  u 


1.2  1 
u  +  u  ;  u 


'4  (  >  U  g  ]  £  )  5 


u(  ; u ^ , Q  ) 


It  then  follows  by  Thns .  3.2,  3.3,  (3.6)— (3.9),  (3.16),  an: 

*  r*  ^ 
the  triangle  inequality  that  u  €  H~  ’  r  (T...  ,T)  and 


(3.1") 


<  ,i,.1m(0,T)  ■■  2.,(7*,T) 

r  +  2,(a  +  2)/2  ~  '  a  1  i  + 2  ,  (  a  +  2 )  /  2  'U  ''  a+  2  ,  (  a  +  2  )  /  2 


<  C{;iuji|a+1  +  'i  ^  /  2  +  T-I(a  +  1)/2 


!u0 !3+l} 


by  (3.14)  and  (3.15)  <  C 


u  i  +  1  c  ;  ^  -  if  3  =  ^ 

UQ '  a  +  1  *  a, a/2’  **  “» 


(l+T,-(a+1)/2)M?R(0,f)  +T.,.-(C+1)/2 


11 U0 11  p+i 

+  ilf : 

;  (  0  ,  T  )  . 

:  a  ,  a/ 2  5 

if  3  <  a_, 

by  (3. 6)- (3.  9) 

S  Cc(p,T*,T 

»u)  . 

It  remains  therefore 

to  construct 

1 

U0 

2 

and  Uq 

satisfying 

( 3 . 11 ) - ( 3 . 15 ) . 

(  Jf  P  =  a, 
(3.18)  / 

set  UJ  *  u0 

and 

u02  =  0, 

and 

if  P  <  a, 

set  uj  = 

and 

u0  =  u0 

-  > 

a 

where  with  i  =  [a/2 

+  1/4]  >  0, 

z . ,  and 

i 

the  sequence 

{z.}.  .  .  is  defined  by: 

]  ]s-l>3- 


z  ^  =  0  and  for  j  =  0,i  z ^  solves 
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where  the  properties  1(6^)  and  IICS^)  are  defined  as 
follows . 

We  say  that  u  €  Ha,^(T.v,T)  satisfies  1(5-^)  if 
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(3.23)  for  each  t  €  [T...  ,T],  ilu  (t,*)i!n  >  0  and 

••  xx  0 


Ca(p,TftlT,u) 


t  €  [ T , T  ]  ,  Uxx  1  ;i|0 


5  V 


We  say  that  u  <  Ha  ’ ^  (T...  ,T)  satisfies  II(5„)  if 

•*  2 


(3.24) 


for  each  t  €  [T.,.  ,T]  3  an  open  interval  I. 

**  t 


such  that 


(3.25)  g.  ( t  ,u) 


in-  <■  <t,x)  i  '  I  I-t- 

x€I  X  T 


}  >  0  and 


(3.26) 


t€[T* ,T]lgL(t,u) 


gLI(t,u)' 


2  f  2  j  where 


(3.27)  g  (t,u) 


sup{|u  (t,x)|  +  I u  ( t ,x) 1 } 
X€I 


Remarks .  From  Thms .  2.2,  3.4,  and  the  Sobolev  Embedding 
Theorem  we  have  that  ( 3 . 23 ) - ( 3 . 27 )  are  well-defined,  since 
for  a  >  1,  u  €  tfa,?(T*,T)  =>  u  €  C°  ( [T,,.  ,T]  ;H2  )  ,  and  for 
a  >  3,  u  €  Ha,1B(TA,T)  -  u  €  C°(  [T*  ,Tl  x  I)  .  Membership 
in  these  subclasses  does  not  depend  on  absolute  magnitude. 


For  u  €  },  it  follows  that  C  •  u  €  } 


61  61,62 


61  61’62 


for  any  nonzero  constant  C  and  the  same  5^{6^  and  Sj}* 
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Let  u(t)  = 
the  solutions  of 


4 .  A  Priori  Estimates 
u(t;uQ,f)  and  U(t)  =  UA(t;U^,PAf)  denote 
eqs.  (2.14)  and  (2.17),  respectively.  Define 


(4.1) 


e(  t ) 

=  u  (  t  )  - 

U(t)  , 

p(  t) 

=  P1u(t> 

-  u(t) , 

9(t) 

=  '■'  P^u(t) 

-  U(t) . 

In  order  to  analyze  the  a  posteriori  error  estimator  in 
t^.me-indepeiWent  space  meshes,  we  require  esti¬ 
mates  of  the  form 


Section  5  to 


1 


(4.2) 


V 


I  9  ( t )  l| r  =  0(h  V 

L 


( e t ( t )  ||Q  =  0(h  n 


as  h 


0;  >  1 ,  ct2  >  0 , 


with  similar  estimates  needed  with  >  1  for  time-dependent 
meshes.  Estimates  of  this  type  for  Galerkin  approximations  of 
time-dependent  equations  have  been  obtained  by  many  authors 
(cf.  [9],  [10],  [17],  [18]),  usually  as  a  means  of  bounding 
some  norm  of  e  for  j  >  0 ,  and  often  where  the  regularity 

of  the  exact  solution  u  is  given  implicitly  by  norms  of  u 
appearing  in  the  estimates.  To  clarify  the  connection  with 
section  5,  throughout  this  section  we  shall  assume  for 
simplicity  that 


(4.3) 


P  ie  a  k(\)  -  regular  mesh  family, 
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and  that  our  choice  of  discrete  initial  data  is  fixed  as 

Kv  1  f  p>  if  uo  !  ®<L)’ 

(4.4)  =  / 

4  e  e,  if  u0  t  »<L>. 

In  a  manner  similar  to  that  used  in  [9],  other  choices 
for  Uq  can  be  handled,  but  we  require  (4.2)  to  hold  uni¬ 
formly  as  t  -*  0  for  our  smooth  data  estimates  valid  on 
[0,T).  Note  that  for  uQ  €  D(L),  U^(0)  =  P^u^CO)  by  (2.9). 

The  main  result  of  this  section  is 

Theorem  4.1.  Let  a,p  satisfy  (3.1)  and  a  <  4. 

Assume  (4.3)  and  (4.4),  let  0  S  T.v  <  T,  and  u(;uQ,f)  € 
Ha,^(Tj.,T).  Then  3  a  positive  constant  C  (depending  on 
the  functions  a  and  b,  constants  X  and  T,  but  not 


on  T.,c 

’  V  f’ 

or 

u)  such  that  for  each  A 

6  P 

(4.5) 

l|0(t)HE 

< 

C  h  X(A)C  (p,T,.{,T,u)> 

(4.6) 

l|et(  t)  II 0 

< 

CT  /  1  t 

C  h  2  (  A )  C  (P,T*,T,U  )J 

a  * 

€  CTft,T), 

where 

(4.7) 

G1 

= 

min(a,2)  , 

f a-K  ,  if  a  s  2 

(4.8) 

= 

L 

1 

1  max( 2-k , a- 2 ) ,  if  2  <  a  < 

4. 
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1  2 

Proof .  For  each  a  and  (3,  let  be  the 

a ,  [3  -dependent  decomposition  of  given  in  (3.18).  Then 

12  3 

u  =  u  +  u  +  u  solves  eqs .  (2.14),  where 

u1  =  u( ;Ug ,f )  , 

(4.9)  (  ii2  =  u(  ;Uq,0)  , 

u3  =  0 . 

V. 

Similarly,  we  decompose  the  solution  U  of  eqs.  (2.17)  as 
U  =  U1  +  U2  +  U3,  where 

V  =  U(iuJ.P0f);  uj 

(4.10)  /  U2  =  U(-,U2,0);  U2 

U3  =  U( ;U3 , 0 )  ;  U3 

We  then  have 
3 

e(t)  =  l  em( t) ;  em(t)  =  um(t)  -  Um(t),  m  =  1,2,3, 
m=l 

p(t)  =  l  pm( t);  pm(t)  =  P1um(t)  -  um(t), 
m=l  1 

0(t)  =  l  8m(t);  9m(t)  =  Pxum(t)  -  Um(t). 

m=l 

Using  the  t-independence  of  the  operators  L,  and 


subtracting  eqs.  (2.16)  from  eqs.  (2.15)  we  get 
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(4.14) 

(4.15) 

(4.16) 

(4.17) 


and  for 

(4.18) 

(4.19) 


£  t  <  T  and  P  5  a 


II  -  n  |u  na  +  2,(a  +  2)/2, 


I  e^"  ( t )  II  <  C  h  2iiu1MC0’T> 

1  0  _  n  ,|U  a+  2  ,  ( a+  2 )  /  2 . 


It  follows  from  (3.18),  (3.20),  (4.4),  (4.7),  (4.8),  (4.10), 
Lemma  2.1,  and  the  triangle  inequality  that  for  p  <  1  <  a 
(i.e.  uQ  |  V(L) ) 

(4.20)  !|Uq|I0  =  ll(p0-p1)tail0  £  C  h2lltaita  +  1  s  C  •  min(h  \h  2) 

M^R(0,f) , 

and  if  1  5  p  <  a  (i.e.  Uq  €  P(L)) 

(^.2i>  iluX  =  »«Vpo)uoMP0-p1)+«"o  5  c  «2<iiu0;!2+ il+all2) 

^9  pn 

<  C  •  min( h  ,h  )  *  {!lu0llp  +  1  +  ArK(0,f)  }. 

Also,  since  i  Ha  + 2  ’  ^ a+ 2^(  G  ,T)  for  p  s  a,  we  have  by 


Thm.  3.2  that 


(4.22) 


flu1: 


( 0  ,T ) 

a+2,U+2)/ 


(  0  ,T ) 
s,"/  2 


}• 


Using  ( 4 . 20 ) - ( 4 . 2 2  )  ,  the  estimates  (2.14)  for 
2 

and  Ug,  and  the  definition  (3.9)  or  in  ( 4 . 14 ) - ( 4 . 19 ) , 

we  see  that  the  demonstration  of  ( 4 . 14 ) - ( 4 . 19 )  will  complete 
the  proof  of  Thm.  4.1. 

We  shall  make  use  of  the  following  simple  lemma. 


Lemma  4.1.  Let  0  2  y^  2  y^  and  0  <  Tft  2  t. 

Then  t  1  2  1  +  T*  2. 

We  first  demonstrate  (4.14)  and  (4.15).  By  (4.9), 
(4.10),  eqs.  (4.11),  and  Lemmas  2.2  and  4 . 1  we  have  that  for 
a  >  1  and  t  2  Tft  >  0 

(4.23)  II  03  ( t )  ltE  =  ||U3(t)|lE  2  C  t"1/2||U3|l0  2  C(l+T;(a+1)/2) 

Kv 

(4.24)  !|e3  ( t)  II  Q  =  ||  ~  U3(t)ll0  <  C  t” 1 II U  3 II  Q  2  C  ( 1+Tl (  a  +  1 }  7  2  ) 

IIUqII  o  5 

which  are  (4.14)  and  (4.15). 

2  2  2 
Because  e  satisfies  eqs.  (4.12)  and  =  PqUq, 

spectral  representations  may  be  applied  to  obtain  (4.16) 

and  (4.17).  A  simple  interpolation  argument  leads  to  the 

following  slight  variation  of  [9,  Thm.  3.1]. 
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Theorem  4.2.  There  exists  a  positive  constant  C  such 
that  for  ]20,0sss2,  and  t  >  T...  >  0 


(4.25) 


0  2  .  ,  v  ,|  r  v,s  ^  s  +  2  j  )  /  2  ,i  2 1 

'at Je  (t)l!n  s  c  h  T*  :|urr 


0  ■'  0 


It  is  easily  verified  in  (4.8)  that 


{ 


a  -  1  €  (0,2),  if  a  €  (1,3), 


2, 


if  a  €  (3,4]. 


The  application  of  (4.25)  with  j  =  1  and  s  =  a-1 
or  2  yields 

(4.26)  II  e^(  t )  i|  Q  <  C  ha_1  T“(  a+ 1}  7  2  ||Uq  II  Q  ;  if  1  <  a  <  3,  and 


(4.27)  ||  e 2  ( t )  ’!  0  <  C  h2  T~2  ||u2 1|  0  ; 


if  3  <  a  s  4, 


The  estimate  (4.17)  then  follows  by  applying  Lemma  4.1. 

2  2  2 

To  show  (4.16)  we  set  p  =  9  =  P  +  e  in  eq.  (4.12) 

and  get  for  t  2  Tft  >  0  and  any  s  €  [0,2],  s  i  1/2 


||  0  2  ( t )  ||p 


<  e2(t),  p2(t)  +  e 2 ( t ) > | 


by  the  triangle  inequality,  (4.1),  and  Lemma  2.1 


<  ||e2(  t)  Il0  {C  hSl|u2(t)  || s  +  l|e2( t)  ll0  } 


by  Lemma  3.1  and  Thm.  4.2 


n  i.2s  .-(S  +  1)|,  2  .j  2 
5  C  h  t  II  u0  II  O’ 


and  hence 
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(4.28) 


9  ( t )  !i 


C  h“ 


T* 


-(s+D/2,  2, 


!iu 


O'O' 


Setting  s  =  =  min(a,2)  for  1  <  a  <  4  and  using 

Lemma  4 . 1  we  get  (4.16). 

We  proceed  to  demonstrate  (4.18)  and  (4.19)  directly 
via  energy  estimates.  These  results  can  also  be  obtained 
for  some,  but  not  all  a.  €  (1,4]  satisfying  (3.1)  by 
applying  Thms .  2.1,  2.2,  and  results  given  in  [17]  in 
arguments  similar  to  (  4 . 25 ) - ( 4 . 28  )  . 

Setting  4)  =  20^(s)  in  eqs .  (4.13)  for  j  =  0  and 
0  <  s  s  t  <  T ,  we  get 

(4.29)  2  i!  01  (  s )  \\2n  +  2B(91(s),  61(s))  =  2<p1(s)  ,01(s)>. 

so  s  s  s 


Using  the  t-independence  of  the  symmetric  form  B  and  the 
fact  that  9(0)  =  0,  integration  of  eq.  (4.29)  on  (0,t) 
followed  by  the  application  of  the  Cauchy-Schwarz  inequality 
yields 

(4.30)  !|01(t)||^  <  ~  f  ||(I-P,  )u^(s)||2  ds, 

E  2  J0  1  5  0 


from  which  we  obtain 


(4.31)  ||01(t)ilE  £  C II  (  I-P^  )u^lj  p0QT  }  . 


For  1  <  a  <  2,  u1  €  Ha+2 ’ ( a+2 } 1  2( 0 ,T)  0  H°((0,T);hJ) 


and  so  by  Lemma  2 . 3  and  (4.7) 


J 


r 
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(4.32) 


II  0 1 C  t )  li  E 


C  h  1||u1| 


(  0  ,T ) 

a+  2 , ( a+  2  )  / 2 


5 


which  is  (4.18)  for  1  <  a  s  2.  The  result  (4.18)  for 
2  <  a  <  4  follows  trivially. 

It  remains  to  demonstrate  (4.19).  For  1  <  a  <  4, 
by  (2.9)  and  (4.10)  we  have  that  for  0  £  t  <  T 


(4.33) 

e^(  t ) 

=  e^(t)  +  e^(t),  where 

(4.34) 

e  ( t ) 

=  ( I-P  )u^( t)  ,  and 

(4.35) 

12,  _ 
et  ( t) 

-•  L  91  ( t )  . 

A 

By  Lemma  2.1,  Thm.  2.2,  and  <  >  1,  we  have  that  for 
1  <  a  :£  2 ,  a  i  3/2 


(4.36)  H  ei1(t)  l|n  £  C  ha-<  sup  ||ui(s)|!  ,  £  C 

11  °  S €  [  0  ,T 3  t  3-1 


a-K 


l1||(0’T) 

a  +  2,  ( a+  2  )  /  2  , 


and  similarly  for  2 
Py  (2.11)  and  (4.18) 

(4.37)  l'e^2  ( t)  II  q  =  ii 


<  a  <  4  with  ha  K  replaced  by  h^  K, 
we  also  have  that  for  1  <  a  <  4 


L91(t)lln  £  C  h_K  II  9 1  ( t )  IL  £  C  h  1 
A  0  L 

1„C0,T) 


a+2 , (a+2 )/2 . 


By  the  triangle  inequality  and  the  definition  (4.8)  of  o  , 
we  see  that  to  complete  the  proof  of  (4.19),  it  suffices  to 


show 
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(4.33)  :|ej(t)!!n  <  C  hU  ^  for  2  <  a  s  4 


^■.nVO’T) 

a  +  2,  (a+2)/2 


We  nave  that  e^( t )  =  (I-P^)u^(t)  +  0^(t) 
It  follows  as  in  (4.36)  that 


(4.39) 


( I-P1 ) u^(  t )  i!  Q  s  C  >aa“ 2  !| u1  M (  0  ’ T } 


a+  2 , ( a+  2 ) / 2 ' 


Setting  o  =  20s(s)  in  eqs.  (4.13)  for  j  =  1  and 


0  <  s  S  t  <  T,  we  get 


(4.40)  -p-i| 0 1  ( s )  t 2  +  2!Ie1(s):i2 

ds  s  0  s  £. 


=  2  <o1  (s) ,  91(s)>, 

ss  s 


Integration  of  (4.40)  on  (Q,t)  followed  by  the  use  of  the 
Cauchy- Schwarz  inequality  yields 


(4.41)  >|  6^(  t )  !j  2 


l|e*(  o)  j! 2  +  c 


o  'i(i-pi>“is<s,||S ds' 


Using  0^(0)  =  ( P. -Pn ) u^( 0) ,  we  get 

c  1  U  L 

(4.42)  j|0*(t)|lo  s  c { !i(P1-i)u^(o) j|0  +  !|(i-P0)uJ(0)!!0  + 

I!  ( i _ p  )U1  ||  ^ o ^  i 
1  tt 11 0,0  }' 


Using  Lemmas  2.1,  2.3,  and  Thm.  2.2  we  get  from  (4.42)  that 


(4.43)  ;|e^(t)ilQ 


C  ha-2!|u1i|^°’2()a+2)/2  for  2  <  a  <  4 . 


The  proof  of  (4.33)  and  Thm.  4.1  is  then  completed  by  using 
(4.39),  (4.43),  and  the  triangle  inequality. 
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Remarks .  The  following  consequences  Thm.  4.1 
relevant  to  the  estimates  in  section  5.  Define  the  i 
<  by 

( (o+l)/2;  1  <  2  <  2, 

(4.44)  <(a)  =  <  3/2  ;  2  <  a  5  3, 

l  2  ;  3  <  c  s  4. 

For  fixed  a  €  (1,4]  and  <  €  [l,<(a)>  set 


(4.45) 


M- C  <  ,c) 


<( a) -< 


It  can  be  verified  from  the  result  of  Thm.  4.1  that 
positive  constant  C  such  that  for  each  A  €  F  and 


t  €  CT...T) 
•*> 


1-K 


(4.46)  h  (A){;|0(t)llE+  i!et(t)i!0} 


C  h24(A)C  (p 


(4.47)  h~<(  A)  !l  0(  t)  i|. 


are 

unction 


—  a 
all 


5.  A  Posteriori  Error  Estimates 


Ail  notation  used  in  Section  4  is  adopted  throughout 

this  section.  The  solution  classes  ^  (T... , T )  and 

^  1 

(T...,T)  are  as  defined  in  (3.21)  -  (3.27). 

'  1  ’  0  2  “ 

We  begin  by  listing  for  easy  reference  the  assumptions 
concerning  solution  regularity,  mesh  family  type,  and 
discrete  initial  data  choice  which  will  lead  to  our  results. 


u  €  Ha  5  ^  ( T... , T )  ;  0  £  T.,.  <  T,  a  and  |3  satisfy 

(5.1)  ^  (3.1)  and  a  £  4.  For  A  €  P,  is  chosen  accordin 


to  (4.4). 


(5.2) 


u  €  H  °|  ’  ^  ( T ,  T )  ;  5  >  0.  P  is  a  <(X)-regular 

°1  _  '  1 
ni6sh  family. 


(5.3) 


u  €  (T.,.,T);  a  >  3,  5,  >  0 ,  6-  >  0 


(5.4)  k  £  Cl,<(a));  <(•)  defined  as  in  (4.44). 


Whenever  they  appear,  unless  otherwise  specified 

(5.5)  etc.  denote  positive  generic  constants  which 

depend  on  the  functions  a  and  b,  constants  X  and 
T,  but  not  on  T{l  ,a,(3 , 6^  ,  §2  ,u,U,  or  the  data  u^ 
and  f . 

We  now  define  certain  concepts  needed  in  our  a  posteriori 
analysis  which  are  related  to  those  developed  for  elliptic 


boundary  value  problems  in  [2],  [3],  [4],  and  [5]. 


For  0  S  T*  <  T  let 


(5.6) 

g :  A  -»■  g 

(5.7) 

i :  A  -*•  ® 

We 

define  the 

(5.8) 

n j (t,gA, iA 

and  the 

associated 

(5.9) 

E(t,g\ 

A  6  P. 


h. 

J. 


- II Li  (t)  -  g  (t)!L  T  ; 

( 12a  .  ;  / ^  °’Ij 


t  €  [T,,,T)  ,  1  <  j  <  N(A) , 


(T  nj(t,g*.^,y« 


t  €  [T.,  ,T)  . 


By  (5.7)  , 

(5.10)  LiA(t,x)  =  -a  (x)iA(t,x)  +  b(x)iA(t,x); 

t , X  €  [T*,T)  X  I  1  <  j  <  N(A) , 

and  therefore  (5.8)  and  (5.9)  are  well-defined. 

Let  $>  €  C^([Tft,T);  H^)  and  g  ,<t>  be  as  in  (5.6)  and 
(5.7).  We  say  that  a  quantity  E(*,g  ,4>)  is  an  upper 
estimator  for  ||<f>  -  4>  il^,  on  [T.,e,T)  if  3  positive  constants 


I 
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is 

Ci,C2,h  such  that  for  A  €  P  with  h(A)  £  C^h*'  and  all 
t  €  [T*,T) 

(5.11)  ||*(t)  -  <^A(t)||E  £  C2  E (t,gA,*A), 


where 

and  C2  are  as 

in  (5.5)  and 

h* 

depends  in 

general  on 

the  mesh  family 

P{ i . e . g .  on 

<} 

and  on  the  class 

in  which 

$  resides  {i.e, 

,g.  on  a,p,51 

» 5  2  } 

• 

If  under  the  same  hypotheses  3  a  constant  a  >  0  such 

that 


(5.12)  ||$(t)  -  /(t)||E  =  E(t,gA,*A)(l  +  0(ha))  as  h  -  0, 

where  the  constant  in  the  0-term  is  as  in  (5.5),  then  we 
say  that  E(*,g,<i>)  is  an  asymptotically  exact  estimator 
for  |(<f>  -  <*>|(E  on  [Tft,T). 

Let  u  and  U  be  the  solutions  of  eqs.(2.14)  and 
(2.17).  We  take  as  our  primary  estimator  for  the  error 
|| u  —  U !! E  the  quantity  E(*,f-Ut,U),  defined  according 
to  (5.8)  and  (5.9).  The  principal  results  in  this  section 
are  the  following  two  theorems. 

Theorem  5.1.  Assume  (5.1),  (5.2),  (5.4),  and  let 
p.  =  p,(K,a)  be  as  defined  in  (4.45).  Then  E(‘,f-U^,U)  is 
an  upper  estimator  for  ||u-U||E  on  [Tft,T),  with  h*  =  5^ 


•■•jcai?'--  <■ 


Jd 


r 


~T 
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Theorem  5.2. 

Assume 

(5.1)  -  (5.4)  and 

let  p  be  as 

defined  in  (4.45). 

Then 

E( * , f-U  ,U)  is  an 

asymptotically 

exact  estimator  for 

!|  u— U 

!Ie  on  [Tft  ,T )  ,  with 

h*  =  min( 5 2^// and  a  -  |i. 


E(*,f-U..,U)  is  based  upon  local,  computable  residuals 
and  can  therefore  be  used  to  monitor  and  control  the  space 
discretization  error  in  an  adaptive  FEMOL  procedure.  While 
h*  may  be  quite  small,  its  existence  indicates  that 
E(*,f  -U^,U)  will  perform  well  as  an  estimator  for  a  given 
class  of  data,  and  not  be  overly  sensitive  to  changes  in  solution 
behavior  during  the  time  evolution  of  a  problem. 

The  proof  of  Thms .  5.1,  5.2  is  based  upon  interpreting 
u^Ct)  as  given  data,  eqs.  (2.14)  as  a  continuous,  one- 
parameter  (t)  family  of  elliptic  boundary  value  problems: 

For  each  t  €  [Tft,T),  find  u(t)  *  P(L)  satisfying 

(5.13)  Lu(t)  =  f(t)-u^.(t), 

and  employing  the  error  I!  p  ( t )  |J and  estimator  E( t , f  -  u^ , P^u) 
associated  with  (5.13).  E(  •  , f  -  u^ , P^u)  is  defined  according 
to  (5.8)  and  (5.9)  and  is  not  computable,  and  ||p|lE  is  not 
of  primary  interest,  but  up  to  higher  order  terms  in  h  we 
have  E(  *  ,f  -  ut,P1u)  ~  E(*,f-Ut>U)  and  ||p|lE~  !)  e  ||  E . 

Specifically,  the  proofs  of  Thms.  5.1  and  5.2  follow 
■J  from  the  application  of  Thm.  4.1  and  the  demonstration  of 

the  following  sequence  of  lemmas. 
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Lemma  5.1.  There  exists  a  positive  constant  C  such 
that  for  each  A  €  P  and  all  t  €  [TA,T) 

lECt.f-U^U*)  -  E(t,f-ut,P*u)|  <  Ch(A){||0(t)!lE  +  II  e  t  ( t )  It  Q}  • 

Lemma  5.2.  Assume  (5.1).  Then  E (• ,f-u  ,P. u)  is  an 
upper  estimator  for  IIp'I^  on  tTsV,T),  with  h*  =  1. 

Lemma  5.3.  Assume  (5.1),  (5.2),  and  the  result  of 
Lemma  5.2.  Then  3  a  positive  constant  C  such  that  for 
each  A  €  P  and  all  t  €  [T* ,T) 

£ 

E(t,f-ut,pju)  >  CA,-1H  Cg(p  ,  T  j,. ,  T ,  u )  . 

Lemma  5.4.  Assume  (5.1)  -  (5.4).  Then  E ( • , f-u^ ,P^u) 
is  an  asymptotically  exact  estimator  for  1!  p  !l  E  on  [T*,T), 
with  h*  =62^  M-  defined  as  in  (4.45),  and  a  =  ijl. 

Before  proving  Lemmas  5.1  -  5.4,  we  show  that  the 
result  of  Thm.  5.1  follows  from  the  validity  of  hypotheses 
(5.1),  (5.2),  (5.4),  and  the  results  of  Lemmas  5.1,  5.2, 

5.3,  and,  if  the  result  of  Lemma  5.4  also  holds,  then  we  have 
the  conclusion  of  Thm.  5.2. 

To  simplify  notation  we  suppress  t,  and  also  write 
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(5.14) 


f'E1  --  E(  t,f-U*,UA), 

<  E ^  =  E(t,f  -  ut,pju), 

=  Ca(0,T*,T,u). 


By  Lemma  5.3  we  may  write 


(5.15) 


1  - 


El"E2 


<  ^  <  ! 
2 


which  after  application  of  Lemma  5.1  and,  again.  Lemma  6.3 
becomes 

C6  h1“IC 

(5.16)  1 - -  £ !|e||E+||et !!q} 

a 


C61h 

— C 


1-K 


a 


(ll9|lE+lletll0>. 


Under  assumptions  (5.1),  (5.2),  and  (5.4),  we  may  apply  the 
consequences  (4.44)  -  (4.47)  of  Thm.  4.1  to  (5.16)  and  get 

(5.17)  1  -  <  ^  <  1  +  CS^h2^. 

Now,  by  the  definition  of  :  Hq  -►  S 


hp«e 


inf  ||u-tl|„  < 
ftS 


9 


(5.18) 
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and  so  by  applying  the  triangle  inequality  to  e  =  - D  +  0 
we  have 


(5.19) 

IIP  II 

E  5  ^ 

< 

E  “ 

^  IIP»E  + 

I|9|Ie. 

By  Lemma  5 . 3 

,  (5 

.19)  becomes 

(5.20)  E2 

IIpIIe 

< 

II e  II E  = 

£  J 

fllpllE  ^ 

CS1h“< 

«ellE], 

E2 

2  1 

L 

C 

a 

and  using  (4 

.44) 

- 

(4.47) 

as 

in  (5.17) 

(5.21)  E2 

1!  ollE 

E2 

< 

II  e  II E  2 

fHp|lE 

L  E2 

csih2  j 

I- 

Setting 

h* 

61  ’ 

we 

see  in 

(5.17  ) 

that  3 

positive  constant  such  that  for  h  £  C^h* , 

(5.22)  E2  =  E-L  (1  +  O(h^) )  as  h^O, 

with  constants  in  O(h^)  as  in  (5.5),  and  therefore  from 
(5.21)  we  get 

l|PllE  u 

(5.23)  ||  e  ||  £  =  E1~ —  (l  +  0(h^))  as  h^O. 

The  demonstration  of  Thms.  5.1,  5.2  now  follows  easily  from 
the  results  of  Lemmas  5.2,  5.4,  respectively. 

We  proceed  to  prove  Lemmas  5.1  -  5.4. 
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Proof  of  Lemma  5.1.  To  further  simplify  notation  we 
write  uA  =  PAu  and  for  1  £  j  2  N(  A) ,  q..(U)  =  ( t , f-UA , UA ) 

and  r^(uA)  =  q ^ ( t , f-u^ ,uA ) .  By  definitions  (5.8),  (5.9), 
and  the  triangle  inequality  we  have 


(5.24) 


where 


\E2-E22\  <  I  |r^(U)-T^(uA)| 

*  j=l  J  J 

N  h2. 

s  rrfr  iQij‘Q2:'  '  iQxjtQ2j4 


Ql.  =  HLU  +  iq-f#^ 

> 

<32j  =  ||LuS+ut-f||0iI^ 


1  2  j  2  N. 


Using  again  the  triangle  inequality,  the  N-dimensional 
Cauchy-Schwarz  inequality,  and  (5.8)  we  get  from  (5.24)  that 


9  ?  f  N  h^  -  9  hl/2 

(5.25)  |e’-E*|  <  2^  5^7 

•  (  I < u)+-n? < ua> > 

v.  j=l  3  J 

It  then  follows  from  (2.1),  (2.2),  (5.9),  and  (5.10)  that 
3  a  positive  constant  C  such  that 


i 
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(5.26)  |E2-E2|  -  Ch(A){||9|lE+||etH0}{E1  +  E2}, 

from  which  the  desired  result  follows. 

Proof  of  Lemma  5.2.  By  Thms .  2.2,  2.3,  and  3.1,  the 
assumption  (5.1)  guarantees  that  u  €  C^([T.V,T);  P(L)). 
Interpreting  eqs .  (2.14)  as  in  (5.13)  and  applying  [2,  Thm. 
(stated  for  u  €  C°(I)  but  valid  in  fact  for  u  €  H°  ; 

XX  '  XX 

cf.  also  [4]),  we  have  the  desired  result. 


Proof  of  Lemma  5.3.  Letting  again 
that  for  t,x  €  [T,..  ,T)  x  I..;  1  2  j  2N 


uA  =  PAu, 


we  note 


(5.27)  LuA(t,x)  -  (f(t,x)  —  ut(t,x)) 

=  (u.-f-a  uA+buA)(t,x)  =  (au  -a  p  +bp)(t,x). 

t  X  X  XX  X  X 

It  then  follows  from  (5.8)  that  3  a  positive  constant  C 
such  that 


(5.28)  ri?(t,f  -  ut,uA)  >  Ch?(||uxx(t,-)||^Ip  -  llp(t)||jsI^}i 


1  <  j  2  N. 


Summing  for  j  =  1,N  yields 

(5.29)  E2(t,f  -  ut,uA)  +  Ch2(A)!|p(t)||2  >  Ch2(A)||uxx(t,*)!!2, 


.1] 


which  by  the  <-regularity  of  assumption  (5.2)  and  the  result 
of  Lemma  5.2  becomes 


(5.30)  E2(t,f-u  ,uA)  >  ( - -  \  h2<(A)||u  (t,.)||2 

\ 1+Ch  ( A)  /  xx  0 

C 

Since  by  (5.2),  Hu^Ct ,  • )  ||  Q  ^  Vt€  CT,.{,T),  we  have 

from  (5.30)  the  result  of  Lemma  5.3. 

We  omit  the  proof  of  Lemma  5.4,  as  it  is  completely 
analogous  to  that  in  [2,  Thm.  7.3],  where  a  sequence  of 
six  asymptotic  equalities  related  to  a  decomposition  of  the 
indicators  led  to  the  desired  result  (cf.  [8]  for  further 
details).  Roughly  speaking,  the  properties  (3.24)  -  (3.27) 
assumed  for  asymptotic  exactness  require  that  the  solution 
is  smooth  and  not  relatively  flat  at  any  time  t  €  [Tft,T). 
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6.  A  Posteriori  Error  Estimates  Revisited 


The  purpose  in  using  a  posteriori,  rather  than  a  priori 
estimates  to  bound  the  error  due  to  the  numerical  solution 
of  differential  equations  is  to  obtain  more  realistic  bounds. 
We  expect  that  restrictive  assumptions  on  admissible  solution 
classes  and  mesh  partitions  need  to  be  imposed  in  order  to 
guarantee  that  an  estimator  is  asymptotically  exact.  However, 
it  is  desirable  that  an  estimator  is  reliable,  i.e.  an  upper 
estimator,  under  far  less  restrictive  and  more  easily 
verifiable  conditions. 

We  shall  show  that  if  the  estimator  E( • , f  -  ,U) 

introduced  in  Section  5  is  modified  through  the  addition  of 
a  term,  which  takes  more  fully  into  account  the  t-dependence 
of  the  problem  and  is  of  higher  order  in  h,  then  we  obtain 
a  reliable  estimator  for  a  class  of  problems  occurring  often 
in  practice,  with  no  restrictions  on  the  mesh  family  or  mesh 
size . 


We  assume  that  the  solutions  u(  ;ug,f)  and 
U(  ;U0,PQf )  of  eqs.  (2.15),  (2.16)  are  such  that 

(6.1)  f,ft  €  C°([0,T);H°), 


(6.2) 


u0  €  Z7(L),  and 


ni 

uo  *  Vo- 


(6.3) 


for  A  €  P. 
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Note  that  from  (6.1)  and  (6.2)  we  have  that 

u  €  C° ( [ 0 ,T ) ;  P(L) )  and  u  €  C°([0,T);  H°). 

2  II *  lip 

Define  C  =  arr  (a  as  in  (2.1))  <  X  =  inf  - . 

4>*  0 

For  t  €  [0,T)  and  A  €  V  we  define  the  estimator 


(6.4) 


E  ( t ,  f  -  UA,UA) 


\  E^ ( t , A )  ,  where 


1  =  1 


(6.5) 


E-^  (t ,  A) 


E(t,f  -  UA,UA) 


ft  * 

(6.6)  E_  ( t ,  A )  =  n  (  A )  {  [  e~C ( t_s 5 E2 ( s , f o -UA  ,UA)ds}1/2, 

*  J  n  S  S  S  S 


(6.7) 


E3(t, A) 


e~C  t/2  •min(l,h(A)/t1/2)  • E( 0 ,UA( 0 )  +  LuQ-f ( 0 ) ,  0  ) 

,  if  t  >  0, 

E(0,UA(0)+LuQ-f(0) ,0) ,  if  t  =  0. 


By  (6.1),  (6.2),  and  the  fact  that  UA  €  C°°((0,T);  S(a))  we 
see  that  E(t,f-UA,UA)  is  well-defined  and  computable, 
given  UA.  Note  that 
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(6.8) 

e2( t,  A) 


N(A)  h. 


jil  12ai 


e’C(t_S)ilA(Us  +  LljA-f)(s’,)|!0  I  ds 


2 

0(h  (A))  as  h(A)  -*■  0,  provided  that  u  is 
sufficiently  smooth  , 


and  that  by  (2.9),  (6.1)  -  (6.3),  and  Lemma  2.1 

(6.9)  E(0,uJ(0)+LuQ-f(0),0) 

N(A)  Ik  1/2 

=  uh  ll(I-Po>CLu0-f(0))K^I_} 

£  Ch (  A)  {  ||Uq  ||  2  +  ||  f  (  0  )  ||  g  }  . 

Suppose  that  for  e  >  0  arbitrarily  small,  we  either 
assume 


(6.10)  in  addition  to  (6.1),  (6.2)  that  6  H2  +  e  and 

f ( 0 )  €  He ,  or  that 

(6.11)  t  >  T*  >  0  and  h(A>  £  Ti/2(1"e)‘ 

Then  from  (6.5),  (6.7)  -  (6.9)  we  see  that 

(6.12)  E(t,f-UA,UA)  =  E(  t ,  f-U A  ,UA)  +  0(h1  +  e  (  A) )  as  h(A)  0. 

Since  the  exponential  decay  rates  in  E2(t,A)  and 
E^(t,A)  are  overly  pessimistic  for  any  given  problem,  and 
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since  E2(t,A)  is  related  to  the  t-dependence  of  the  local 

indicators  (h  •  }  •  ,  ,  .  s  ,  a  practical  strategy  for  the 

]  3  =  1 , i.  ; 

implementation  of  E(t,f-U  ,U)  would  be  to  use  only 

_  • 

E(t,f-Ut,U),  andto  monitor  some  measure  of 

on  a  time  interval  [t-CghY,t],  for  some  positive  constants 
Cq  and  y-  We  give  no  details  here. 

We  shall  show  the  following  result  . 

Theorem  6.1.  Assume  (6.1)  -  (6.3).  Then  E(*,f-U  ,U) 
is  an  upper  estimator  for  Hu-Uilj,  on  [0,T),  with  h*  =  1. 

Proof.  The  proof  is  based  on  a  different  decomposition 


of  the  error 

e  =  u 

-  U  than  that  used  in  Section  5.  We 

write 

(6.13) 

e(t) 

=  e1(t)  +  e2(t),  where 

(6.14) 

e^  ( t ) 

=  vA(t)  -  UA(t), 

(6.15) 

e2(t) 

=  u(t)  =  vA(t)  , 

and  vA  € 

(  [  0  ,  T )  ; 

P(L))  is  defined  by 

(6.16) 

LvA(t) 

=  (f-UA)(t);  t  6  [ 0 , T ) . 

By  (6.1)  and 

the  t- 

independence  of  the  operator  L,  we 

have  that 

vA  €  C° ( [ 0 ,T)  ;  P(L)) 

and 


£0 


(6.17)  LvJ(t)  =  (f . -uj . )(t)  ;  t  €  [0,T). 

L  L  u  u 

From  eqs ,  (2.15),  (2.16),  and  the  definition  of  P^:  Hg  -*  S(  l 1) 
it  follows  that 

(6.18)  P *  v ( t ,  • )  =  ~-UA(t,-);  t  6  [ 0 , T)  ,  j  =  0  or  1. 

Subtraction  of  eqs.  (2.16)  from  eqs.  (2.15)  shows  that 


(6.19)  e  (t) 


21  22 

e  ( t )  +  e  ( t ) ,  where 


<e^ll(s),i>  +  3(e21(s),i>)  =  -<e1  (  s  )  ,*> ;  Vs  ,<<>€(  0  ,T )  *  hJ  , 


(6.20) 


21 

e^(0)  =  0, 


r<el2(s)  ,*>  +  B(e22(s )  ,4>)  =  0;  V  s,«  6  (0,T)  *  hJ  , 
s  u 


(6.21)  { 


e2  2 ( 0  )  =  e2(0). 


We  show  that  there  exists  a  positive  constant  C  such  that 
for  all  t  6  [0,T)  and  each  A  €  P 

(6.22)  He1  (t )  f[E  ^  C  EL(t,A)  , 

(6.23)  ||e21(t)||E  <  C  E2(t,A) ,  and 

(6.24)  lie2  2  ( t )  ||E  <  C  Eg  ( t ,  A )  , 

which  will  complete  the  proof  of  Thm.  6.1.  By  (6.16)  and 

(6.18),  the  estimate  (6.22)  follows  from  the  result  in 
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[2,  Thm.  7.1],  as  in  the  proof  of  Lemma  5.2. 

From  the  t-independence  of  the  operator  L  and 

^1S/^  21 

eqs.  (6.20),  we  get  that  the  function  w(s)  =  e  e  (s) 

satisfies 


(5.25) 

X  \s/2 

^<w  ( s ) ,  i>  +  B(w(s),V>)  -  —  <w(s),<5>  =  -e  <e^(s),<5>; 


Vs €  ( 0  ,T)  x  Hj, 


w(  0 )  = 


Setting  <t>  =  2wg(s)  in  ea.  (6.25)  yields 


(6.26)  2!lws(s)U*  +  ^  f *|w(s)!l|  -  ||w(s))!g} 

X,  s/2 


-2e  ^  <e^(s),w  (s)>. 

s  ’  s 


Integrating  (6.26)  on  (0,t)  for  t  €  (0,T),  applying 
the  Cauchy-Schwarz  inequality,  and  using  w(0)  =  0  yields 


(6.27)  ||w(t)!|£  -  —  (Iw(t )  ||  q  <  \ 


t  X  s  i  2 

e  1  ||e^( s )  || fds  , 
0  s  u 


which  by  the  definitions  of  w,X^,  and  C  becomes 


(6.28) 


|e21(t)||E  <  { 


■e-c(t-s)||el(s)||2ds}1/2 


By  (6.18),  B(e1(s),t)  =  0  Vs,t  €  (0,t)  x  S(A),  and  so 


by  a  standard  duality  argument  we  have  that 
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(5.29)  ||  e  21  ( t )  ||  £  <  Ch(  A)  { 


'  -C(t-s),,  1,  .,,2,  ,1/2 
e  tl e s  (  s  )  ||  Eds  } 


Using  the  result  in  [2,  Thm.  7.1],  as  in  the  demonstration 
of  (6.22),  we  get  (6.23)  from  (6.29). 

\ls/2  22 

Defining  now  z(s)  =  e  e  “"(s),  from  eqs.  (6.21)  we 

get 


(6.30) 


<zs(s )  ,<t>>  +  B(z(s),<j>)  -  -j-  <z(s),«>  =  0; 


Vs  ,<*  €  (  0  ,T)  x  Hj  , 


z ( 0 )  =  e  (0)  . 


Setting  4>  =  2z  (s)  in  (6.30)  and  using  the  same 
s 

arguments  which  led  to  (6.28)  yields 


(6.31) 


|e22(t)|!E  <  21/2e'Ct/2||e2(0)||E. 


Consecutively  setting  <f>  =  z(s)  and  2z  (s)  in  (6.30) 

S 

and  integrating  on  (0,t)  yields 

(6.32)  2t{||z(t)||2  -  ||z(t)||2}  <  II  z(  0 )  ||  q  , 


which  by  the  definitions  of  z,  ,  and  C  becomes 

-Ct/2 


(6.33) 


|e22(t)||r,  <  - 


27  2  He  <  0  >  ||  Q  ,  for  t  >  0, 


It  follows  from  (6.3),  (6.14),  and  (6.15)  that 
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(6.34)  B  ( e 2  ( 0 )  ,\|r )  =  0  V\|r  €  S(  A )  , 

which  again  by  a  standard  duality  argument  shows  that 

(6.35)  ||  e  2  (  0  )  ||  Q  <  Ch  (  A  )  |j  e  2  (  0  )  ||  £ , 
and  by  the  result  in  [2,  Thm.  7.1]  that 

(6.36)  ||e2(0)||E  S  CE(  0  ,Le2  (  0 )  -0 , 0 ) 

=  CE(0,U*(0)+LuQ-f (0) ,0) . 

Combining  (6.31),  (6.33),  (6.35),  and  (6.36),  we  have  (6.24), 
thus  completing  the  proof  of  Thm.  6.1. 
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7 .  Computational  Procedures  and  Examples 

In  this  section  we  outline  the  main  features  of  a 
general  FEMOL  program  which  includes  a  posteriori  error 
analysis.  The  program  is  of  a  research  type  for  assessing 
various  aspects  and  performances  of  the  FEMOL  procedure. 

While  it  is  not  presently  for  commercial  use,  it  includes 
many  user-oriented  features  employed  in  available  commercial 
software.  The  experiments  using  the  program  were  directed 
toward  the  evaluation  of  the  efficiency  of  the  approach  and 
the  applicability  of  conclusions  based  upon  the  asymptotic 
analysis  presented  in  Sections  4-6. 

The  discussion  here  is  principally  oriented  to 
selected  examples,  in  which  the  exact  solutions  are  available. 
These  solutions  have  forms  typical  of  those  arising  often 
with  systems  of  linear  and  nonlinear  parabolic  equations 
used  in  various  applications,  and  include  decay  and  sharp 
transitions  in  time,  oscillations,  and  travelling  waves. 

The  model  problem  considered  is  that  given  by  eqs.  (1.1), 
where  the  coefficients  a(x)  =  cosh(4x-2)  and 
b(x)  =  sinh(2x)  are  as  pictured  in  figure  2,  and  the 
functions  u^,  f,  and  g  are  chosen  to  be  the  smooth, 
compatible  data  such  that  the  exact  solution  u  of  eqs.  (1.1) 
is  as  given  in  each  of  the  examples  below.  The  theory  in 
sections  2-6  was  presented  for  the  case  when  the  boundary 
data  g  =  0,  but  holds  more  generally  for  eqs.  (1.1). 
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For  t  >  0,  our  goal  is  estimate  !||e(t,*)||!  with  E(t), 
where  e  =  u  -  U  is  the  space  discretization  error,  the 
computable  a  posteriori  estimator  £(•)  is  as  defined  in 
(1.4)  and  (5.9),  and  the  norm  |(|  *  ||!  is  as  defined  in  (1.2). 
The  theory  in  sections  2-6  was  presented  for  the  estimation 
of  ||e(t,*)[lp,  where  the  norm  l!*|lr  is  defined  in  (2.4), 
but  we  note  that  up  to  higher  order  terms  in  the  space  mesh 
size,  !|e(t ,  •  )||£  ~  III  e(t,  •  )|||  . 

Our  evaluation  of  E(*)  is  based  upon  the  behavior  of 
the  effectivity  ratio  0(*),  defined  in  (1.5),  and  we  shall 


examine  how  it  depends  on  the  relative  error 
defined  in  (1.9). 


E  (  •  ) 
^REL/  ’  ’ 


as 


The  theory  given  in  Section  5  showed  that  for  sufficiently 
(very)  fine  partitions  of  the  interval  (0,1),  l/0(t)  is 

bounded  uniformly  for  all  t  in  a  time  interval  of  interest, 
and  0(t)  -*•  1  as  the  partitions  are  refined.  We  shall  see 
that  0(t)  is  near  1  for  practically  any  partition  of  (0,1) 
in  the  examples,  and  in  fact  appears  to  converge  very  rapidly 
to  1  as  the  partitions  are  refined. 

Before  presenting  our  results,  we  first  discuss  some 
of  the  relevant  aspects  of  the  computational  procedures  used 
in  the  experiments  and  more  general  problems.  For  further 
details  we  refer  to  a  following  paper  [7]  and  [8]. 

All  computations  were  carried  out  in  double  precision 
arithmetic  by  the  research  program  FEM0L1 ,  which  was 
written  and  implemented  on  the  IBM  370  System  in  the  Division 
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of  Computer  Pesearch  and  Technology,  NIH.  FEM0L1  has  the 
capabilities  to  solve  coupled  systems  of  linear  and  nonlinear 
parabolic  PDEs  with  piecewise  linear  elements  in  one  space 
dimension,  with  smooth  or  unsmooth  initial  data,  and  general 
separated  end-point  boundary  conditions.  The  program  has 
four  basic  stages  of  operation,  consisting  of  Cl)  the  processing 
of  control  information,  (2)  the  assembly  of  matrices,  (3)  the 
initialization  of  data,  and  (4)  the  time-integration  of  the 
resulting  ODE  system.  Many  of  the  standard  user-oriented 
features  and  logical  switches  employed  in  these  four  stages 
were  borrowed  from  existing  software,  much  of  which  was 
surveyed  in  [16].  When  discontinuous  changes  in  the  space 
mesh  are  allowed,  decisions  concerning  these  changes  are  made 
adaptively  and  looping  occurs  in  stages  2-4.  We  shall  restrict 
most  of  the  discussion  to  those  features  relevant  to  the 
experiments  conducted  here. 

In  the  control  stage,  flags  are  set  which  govern  the 

number  and  type  of  PDEs  solved,  the  space  mesh  used,  and  whether 

or  not  the  error  is  to  be  estimated,  the  true  error  is  to  be 

computed,  or  adaptive  mesh  construction  is  to  be  employed. 

Here,  the  target  time  points  for  output  are  also  determined. 

In  the  examples  to  be  presented,  100  eaually  spaced  time 

points  {T  }  ,  ,nn  were  distributed  over  the  time  intervals 

m  m=l ,iuu 

of  interest,  at  which  E,  |||  e  ||j  ,  and  ®  were  computed, 

but  other  selections  for  output  are  possible. 
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As  in  the  computation  of  E,  |||  e  |||  ,  and  ®» 

FEM0L1  uses  standard  two-point  Gaussian  quadrature  in  the 
assembly  stage  to  compute  the  mass  and  stiffness  matrices 
and  load  vector  needed  for  the  reduction  of  eqs.  (1.1)  to  an 
implicit  system  of  ODEs.  This  process  of  course  introduces 
the  effects  of  the  numerical  integration  into  the  computations. 
After  minimal  testing,  this  effect  was  accepted  for  our 
experiments  with  time-independent  space  meshes.  In  some  of 
the  computations  performed  with  adaptively  constructed  meshes, 
however,  these  pollution  effects  were  more  significant  and 
an  alternate  procedure  had  to  be  adopted  (cf.  a  following 
paper  [7]). 


With  the  smooth,  compatible  data  in  the  examples  here, 
the  discrete  initial  data  was  determined  in  the  initialization 
stage  as  the  standard  finite  element  approximation  of  the 
solution  of  an  elliptic  boundary  value  problem.  When  the 
initial  data  is  not  smooth,  FEM0L1  chooses  the  linear 
interpolation  as  the  discrete  approximation. 

The  time-integration  of  the  resulting  ODE  system  is 
accomplished  in  a  subroutine  with  a  version  of  Gear’s 
variable-step,  variable-order  backward  differentiation 
formulas  which  we  modified  to  efficiently  handle  matrices 
with  banded  structure.  The  ODE  solver  is  essentially 
equivalent  to  that  in  [12].  To  ensure  that  the  effects  of 
the  time  discretization  did  not  pollute  the  results,  very 
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stringent  error  requirements  were  placed  on  the  time 
discretization  scheme.  The  ODE  solver  used  places  equal 
weights  on  the  N-l  components  {U(t,Xj)}j_^  ^  ^  in 
attempting  to  control  the  relative  time  discretization  error 
per  step.  These  are  not  the  best  criteria,  in  that  each 
component  is  given  equal  weight  and  that  an  error  per  unit 
step  tolerance  probably  has  more  meaning.  While  some  minor 
modifications  in  this  procedure  were  made,  these 
inefficiencies  were  accepted  for  the  sake  of  illustrating 
the  theoretical  results.  We  experimentally  adjusted  the  input 
tolerance  until  the  discrete  time  derivative  U  in  each 

_3 

example  changed  by  less  than  10  throughout  the  time 

intervals  of  interest.  This  required  relative  error  per 

-  7  -8 

step  tolerances  on  the  order  of  10  to  10  .  We  emphasize 

that  these  stringent  requirements  should  not  be  used  in 
practice,  but  were  made  here  to  isolate  the  effects  to  be 
sti died . 

One  of  the  main  purposes  in  studying  a  posteriori  error 
estimates  in  the  FEMOL  solution  of  PDEs  is  the  development 
of  more  efficient  algorithms.  However,  because  of  the 
stringent  requirements  set  above  and  the  mentioned,  and  other 
not  mentioned  inef ficiences  which  were  accepted  for  our 
convenience,  a  quotation  of  CPU  times  used  in  the  examples 
is  not  relevant. 


As  a  final  computational  note,  we  remark  that  all  of 
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the  logical  switches  and  capabilities  of  FEM0L1  can  be  set 
for  a  given  problem  by  supplying  subroutines  for  the  problem 
data  and  modifying  less  than  20  cards  of  input  data. 

Example  1.  For  0  2  t  2  3  the  exact  solution  u  of 
eqs.  (1.1)  is  taken  to  be 

(7.1)  u(t,x)  =  u1(t,x)  +  [~  +  j tanh(lOt-ll)]  •  u2(t,x), 
where 

(7.2)  u1  ( t ,  x )  =  1  -  exp ^ ,  and 

2 

(7.3)  u  (t,x)  =  2  sin(irx)  +  2  cos(2irt)sin(2Trx) . 

For  0  2  t  2  .9  u1  completely  dominates  the  behavior 

of  u  and  u  decays  in  t  (cf.  figure  3(a)).  For 

.9  <  t  <  1.1  u  undergoes  a  sharp  transition,  and  for 

2 

1.1  2  t  2  3  the  oscillatory  character  of  u  is  dominant 
(cf.  figure  3(b)). 

The  FEMOL  was  implemented  in  this  example  for  a 
sequence  of  three  uniform  space  meshes,  using  10,  20,  and 
40  finite  elements.  Figure  4  shows  the  behavior  of  the  errors 
in  time,  and  the  linear  rate  of  convergence  is  apparent. 

For  t  >  1  the  errors  in  this  sequence  oscillate  about 
the  10%,  5%,  and  2.5%  levels.  Figure  5  illustrates  how 
well  the  estimator  £(•)  works.  For  any  of  the  meshes  and 
all  t  €  CO, 3]  we  see  that  the  effectivity  ratio  ©(•) 
satisfies 
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(7.4)  1 0  C  t )  -  1 |  £  .07, 
and  that  with  40  elements 

(7.5)  1 0(t )  -  1 |  £  .  005  ;  Vt  €  [0,3]. 

In  addition  to  these  unexpectedly  good  results,  it 
appears  that 

(7.6)  0(t)  =  1  +  0(N-2)  as  N  h.  vt  £  [0,3], 

a  result  much  better  than  that  which  we  could  prove  in 
Thm.  5.2. 

Example  2.  For  0  £  t  £  .1  the  exact  solution  of 
eqs.  (1.1)  is  taken  to  be  that  defined  in  (1.8),  namely 

(7.7)  u(t,x)  =|+|  tanh  C2p(x-10t)];  p  =  20. 

The  solution  is  a  wave  with  approximate  front  width 
p-"*-  which  moves  in  the  positive  x-direction  at  speed  10, 
and  is  as  pictured  in  figure  1  in  section  1. 

As  in  Example  1,  the  estimator  £(•)  was  tested  by 
implementing  the  FEMOL  for  a  sequence  of  uniform  space 
meshes,  using  here  20,  40,  80,  and  160  finite  elements.  In 
figure  6  we  see  the  oscillatory  temporal  character  of  the 
error  for  20  and  40  elements,  due  to  the  inability  of  the 
mesh  to  resolve  the  wave  front.  The  error  again  appears  to 
converge  linearly  for  all  t,  provided  that  N  ?  40.  Note 
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that  for  N  =  20(40)  the  error  oscillates  about  the  50%  (25%) 
mark.  In  figure  7  we  see  that  even  in  the  case  of  50%  error, 
the  effectivity  ratio  0(*)  satisfies 

(7.8)  j  <  0 ( t )  <2;  Vt  e  [ 0  , .  1  ]  . 

Figure  8  illustrates  that  with  an  average  error  of  about 
25%,  for  almost  all  t  €  [0,.l] 

(7.9)  | © (t ) -1 |  5  .1, 

and  that  when  N  =  160  and  the  error  is  about  6%  for  each 
t,  we  have 

(7.10)  |0(t)-l|  <  .005;  Vt  6  [0,.l]. 

The  quality  of  these  results  does  not  appear  to  be 
limited  to  the  use  of  uniform  space  meshes,  as  experiments 
with  slightly  nonuniform  meshes  in  the  same  examples  has 
indicated.  Also,  we  shall  see  in  a  following  paper  [7]  that 
the  adaptively  constructed  meshes  for  Example  2  can  be  very 
nonuniform,  and  yet  the  analogous  effectivity  ratio  is  close 
to  1 . 

We  finally  remark  that  the  norms  il'lljr  and  III*  III  used 
in  the  analysis  and  computations  were  chosen  because  they  are 
the  natural  norms  arising  from  the  weak  form  of  eqs.  (1.1). 

In  applications  where  there  is  no  natural  energy  associated 
with  the  system,  or  when  some  other  characteristic  is  of 
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primary  interest,  then  it  may  be  desirable  to  estimate 

a  different  norm  of  the  error.  By  altering  the  approximating 
finite  element  subspaces  and  a  posteriori  estimator  in  our 
procedure  appropriately,  such  estimates  can  be  obtained. 
Nonetheless,  the  above  results  illustrate  the  robustness  of 
the  estimator  EC*),  even  when  apparently  out  of  the 
asymptotic  range. 
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Figure  2.  Coefficients  of  Eqs.  (1.1)  vs.  X 


u 

(a) 


U 

(b) 


Figure  3 
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X 


Exact  Solution  u(T,X)  of  Example  1  vs.  X 
for  Various  T. 


Figure  6. 


Relative  Error  EREL(T)  in  Example  2  vs.  T 
N  =  Number  of  Uniform  Finite  Elements. 


Figure  7.  Effectivity  Ratio  0(T)  in  Example  2  vs.  T. 
N  =  Number  of  Uniform  Finite  Elements. 


Effectivity  Ratio  0(T)  in  Example  2 
vs.  T.  N  =  Number  of  Uniform  Finite 
Elements . 
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